
############################################
############################################
# Replication file:                       ##
# Conceptions of National Identity        ##
# and  Ambivalence towards Immigration    ##
# Lindstam et al. (2018)                  ##
############################################
############################################

# data files needed: data2.Rdata, data4.Rdata, dataQ.Rdata, dataS.Rdata

##################################
# Analysis with 2015 YouGov data #
##################################

# Load data
load("data2.Rdata")

#######################
#######################
# Model 1 of Table 3 ##
# + Figure 1         ##
#######################
#######################

# Ordered regression
library(MASS)
library(survey)
library(AER)

# Reults with survey package
design<-svydesign(id=~1, weights=~weight, data=data2)
modelord <- svyolr(as.ordered(ref3)~civ+ethn+attach+pint+treat3+pi_afd+pi_fdp+pi_greens+pi_left+pi_spd+pi_cducsu+age+edulo+eduhi+female+ost, design=design)
summary(modelord)
coeftest(modelord)

###########
# PLOTS   #        
###########

attach(data2)

ethnic <- seq(0,1, length.out=100)
civic <- seq(0,1, length.out=100)

# Scenario 1 (range of ethnic identity)
shortdata <- cbind(mean(civ),ethnic,mean(attach), mean(pint), 0, mean(pi_afd), mean(pi_fdp), mean(pi_greens),
             mean(pi_left), mean(pi_spd), mean(pi_cducsu), mean(age), mean(edulo), mean(eduhi), mean(female), mean(ost))

# Calculate linear predictor etc.
ystar <- shortdata%*%coef(modelord[1:16]) 
m <- 5
tau <- modelord$zeta
probs <- matrix(nrow=length(ystar),ncol=m)

# Predicted Probailities
probs[,1] <- pnorm(tau[1]-ystar) # 1st category
for (j in 2:(m-1)) # categories inbetween
probs[,j] <- pnorm(tau[j]-ystar) - pnorm(tau[j-1]-ystar)
probs[,m] <- 1-pnorm(tau[m-1]-ystar) # last category

# Scenario (range of civic identity)
shortdata2 <- cbind(civic, mean(ethn),mean(attach), mean(pint), 0, mean(pi_afd), mean(pi_fdp), mean(pi_greens),
              mean(pi_left), mean(pi_spd), mean(pi_cducsu), mean(age), mean(edulo), mean(eduhi), mean(female), mean(ost))

# linear predictor 
ystar2 <- shortdata2%*%coef(modelord[1:16]) 
probs2 <- matrix(nrow=length(ystar),ncol=m)

# Predicted Probailities
probs2[,1] <- pnorm(tau[1]-ystar2) # 1st category
for (j in 2:(m-1)) # categories inbetween
probs2[,j] <- pnorm(tau[j]-ystar2) - pnorm(tau[j-1]-ystar2)
probs2[,m] <- 1-pnorm(tau[m-1]-ystar2) # last category

# Prepare for plot with GGplot
library(ggplot2)
library(reshape2)
df <- rbind(probs,probs2)
#df <- as.data.frame(df)
plotdf <- as.data.frame(df)
head(plotdf)
names(plotdf) <- c("1-Agree", "2", "3", "4", "5-Disagree")
identity <- rep(1:2, times=1, each=100)
plotdf <- cbind(melt(plotdf),ethnic,identity)
names(plotdf) <- c("Value", "value", "ethnic","identity")
head(plotdf)

labels <- c("1" = "Ethnocultural", "2" = "Civic")

# Make plot
#pdf("plot1.pdf")
ggplot(plotdf, aes(x=ethnic,y=value, fill=Value,group=Value)) +
  geom_area(position = 'stack', colour="black",alpha=0.5) +  scale_fill_brewer(palette = "Blues")+
  ggtitle("Effect of National Identity Conception", subtitle="On Support for Helping Refugees (1-5)") + 
  labs(x="Identity Conception", y="Predicted probabilities") + facet_wrap(~identity, labeller=labeller(identity=labels)) +  
  theme(aspect.ratio = 1.2)
#dev.off()

############################
############################
# Models 2,3,4 of Table 3 ##
# + Figure 2              ##
############################
############################

# Load oglmx - package for estimating heteroskedastic model (Carroll, 2016).
library("oglmx")
load("data4.Rdata")

# Main model (model 2)
modelhet <- oglmx(ref3~civ+ethn+attach+pint+treat3+pi_afd+pi_fdp+pi_greens+pi_left+pi_spd+pi_cducsu+edulo+eduhi+age+female+ost,
                  ~civ+ethn+civ*ethn+attach+pint+treat3, data=data2, link="logit", weights=weight, constantMEAN=FALSE, constantSD=FALSE)

summary(modelhet)

# Model with additional controls (model 3)
modelhet2 <- oglmx(ref3~civ+ethn+attach+pint+treat3+pi_afd+pi_fdp+pi_greens+pi_left+pi_spd+pi_cducsu+edulo+eduhi+age+female+ost,
                   ~civ*ethn+trad*uni+attach+pint+dist+pi_strength+treat3+pi_afd+pi_fdp+pi_greens+pi_left+pi_spd+pi_cducsu+edulo+eduhi+age+female+ost, data=data4, link="logit", weights=weight, constantMEAN=FALSE, constantSD=FALSE)

summary(modelhet2)

# Model with nominal measures (model 4)
modelhet3 <- oglmx(ref3~civ+ethn+attach+pint+treat3+pi_afd+pi_fdp+pi_greens+pi_left+pi_spd+pi_cducsu+edulo+eduhi+age+female+ost,
                   ~as.factor(dummy)+attach+pint+treat3+pi_afd+pi_fdp+pi_greens+pi_left+pi_spd+pi_cducsu+edulo+eduhi+age+female+ost, data=data2, link="logit", weights=weight, constantMEAN=FALSE, constantSD=FALSE)

summary(modelhet3)

# compare to simple model without modeling the variance 
modelsimple <- oglmx(ref3~civ+ethn+attach+pint+treat3+pi_afd+pi_fdp+pi_greens+pi_left+pi_spd+pi_cducsu+edulo+eduhi+age+female+ost,
                     data=data2, link="logit", constantMEAN=FALSE, weights=weight, constantSD=FALSE)

summary(modelsimple)

modelsimple2 <- oglmx(ref3~civ+ethn+attach+pint+treat3+pi_afd+pi_fdp+pi_greens+pi_left+pi_spd+pi_cducsu+edulo+eduhi+age+female+ost,
                      data=data4, link="logit", constantMEAN=FALSE, weights=weight, constantSD=FALSE)

summary(modelsimple2)

##############
# Model fit  #
##############

library("lmtest")

lrtest(modelsimple,modelhet)
lrtest(modelsimple2,modelhet2)
lrtest(modelsimple,modelhet3)

#############
# Figure 2  #
#############

# Extract coefficients
coef <- coef(modelhet)
coef <- as.matrix(coef)
dim(coef)

# Extract variance-covariance matrix
vcov <- vcov.oglmx(modelhet)
dim(vcov)

# Draw 1000 simulated coefficients 
S <- mvrnorm(1000, coef, vcov)
Ss <- S[,17:22]

##############
# Scenario 1 #
##############

# Conception range
range <- seq(0,1, length.out=1000)

# Set covariates
values1 <- cbind(0.1,range, mean(attach), mean(pint), mean(treat3), 0.1*range)

# Calculate linear predictor  
props1 <- Ss%*%t(values1)

# Means and confidence intervals 
p.mean <- exp(apply(props1,2,mean))
p.qu <- exp(t(apply(props1,2,quantile,prob=c(0.025,0.975))))

##############
# Scenario 2 #
##############

# Set new covariates
values2 <- cbind(0.9,range, mean(attach), mean(pint), mean(treat3), 0.9*range)
                
# Calculate linear predictor 
props2 <- Ss%*%t(values2)

# Mean and confidence intervals 
p.mean2 <- exp(apply(props2,2,mean))
p.qu2 <- exp(t(apply(props2,2,quantile,prob=c(0.025,0.975))))

##############
# Scenario 3 #
##############

# Set covariates 
values3 <- cbind(range,0.1, mean(attach), mean(pint), mean(treat3), range*0.1)

# Calculate linear predictor 
props3 <- Ss%*%t(values3)

# Mean and confidence intervals 
p.mean3 <- exp(apply(props3,2,mean))
p.qu3 <- exp(t(apply(props3,2,quantile,prob=c(0.025,0.975))))

##############
# Scenario 4 #
##############

# Set covariates
values4 <- cbind(range,0.9, mean(attach), mean(pint), mean(treat3), range*0.9)
                
# Calculate linear predictor 
props4 <- Ss%*%t(values4)

# Mean and confidence intervals
p.mean4 <- exp(apply(props4,2,mean))
p.qu4 <- exp(t(apply(props4,2,quantile,prob=c(0.025,0.975))))

#################################
# Prepare the data for the plot #
#################################

# Append the different predictions and CIs
one <- list(p.mean, p.mean2,p.mean3,p.mean4)
one <- unlist(one)
two <- list(p.qu[,1], p.qu2[,1],p.qu3[,1], p.qu4[,1]) 
two <- unlist(two)
three <- list(p.qu[,2], p.qu2[,2],p.qu3[,2], p.qu4[,2])
three <- unlist(three)

# 4 times 0-1 range. 
x <- seq(0,1, length=1000)
s <- rep(x,4)

# 2 times 'low' and 'high' 
Conception <- rep(0:1,2, each=1000)
Conception <- ifelse((Conception==0),"Low", "High")

# Put data together 
identity <- rep(1:2, times=1, each=2000)
df.plot <- cbind.data.frame(one, two, three, s, Conception,identity)

labels <- c("1" = "Ethnocultural", "2" = "Civic")

#pdf("plot2_grey.pdf")
ggplot(df.plot, aes(x = s, y = one, ymax = three, ymin = two, group = Conception)) +
  geom_ribbon(alpha = 0.3, aes(fill = Conception))+ geom_line(aes())  + theme_bw() + scale_fill_manual(name="Other \nConception", values=c("dodgerblue4","skyblue3")) +
  ylab("Standard Deviation of Residual") + xlab("Identity Conception") +
  ggtitle("Effect of National Identity Conception", subtitle="On Response Variability") + facet_wrap(~identity, labeller=labeller(identity=labels)) +  
  theme(aspect.ratio = 1.2)
#dev.off()

############################
############################
# Malleability:           ##
# Table 5 + Figure 4      ##
############################
############################

# ols without controls - main effect 
model1 <- svyglm(ref3~treat3, design=design)
summary(model1)

# ols without controls, heterogeneous effects 
model2 <- svyglm(as.numeric(ref3)~as.factor(dummy)*treat3, design=design)
summary(model2)

# ols with controls, heterogeneous effects 
model3 <- svyglm(as.numeric(ref3)~as.factor(dummy)*treat3+attach+pint+pi_afd+pi_fdp+pi_greens
                 +pi_left+pi_spd+pi_cducsu+age+edulo+eduhi+female 
                 +ost, design=design)
summary(model3)

# ordered logistic, heterogeneous effects 
model4 <- svyolr(as.factor(ref3)~as.factor(dummy)*treat3+attach+pint+pi_afd+pi_fdp+pi_greens
                 +pi_left+pi_spd+pi_cducsu+age+edulo+eduhi+female 
                 +ost, design=design)
summary(model4)
coeftest(model4)

# ordered logistic, heterogeneous effects - reference category: ethno-cultural  
model5 <- svyolr(as.factor(ref3)~as.factor(dummy2)*treat3+attach+pint+pi_afd+pi_fdp+pi_greens
                 +pi_left+pi_spd+pi_cducsu+age+edulo+eduhi+female 
                 +ost, design=design)
summary(model5)
coeftest(model5)

# ordered logistic, heterogeneous effects - reference category: civic 
model6 <- svyolr(as.factor(ref3)~as.factor(dummy3)*treat3+attach+pint+pi_afd+pi_fdp+pi_greens
                 +pi_left+pi_spd+pi_cducsu+age+edulo+eduhi+female 
                 +ost, design=design)
summary(model6)
coeftest(model6)

##################
# Make plots     #
##################

# Extract coefficients
coef <- coef(model3)
coef <- as.matrix(coef)
dim(coef)

# Extract variance-covariance matrix
vcov <- vcov(model3)
dim(vcov)

# Draw 1000 simulated coefficients 
S <- mvrnorm(1000, coef, vcov)

##############
# Scenario 1 #
##############

# Set covariates: Mixed conceptions, no treatment 
covar <-  cbind(1, 0,1, 0, mean(attach), mean(pint), mean(pi_afd), mean(pi_fdp),
                mean(pi_greens), mean(pi_left), mean(pi_spd), mean(pi_cducsu), mean(age),
                mean(edulo), mean(eduhi), mean(female), mean(ost), 0,0)

# Calculate predictons 
props1 <- S%*%t(covar)

# Means and confidence intervals 
p.mean <- apply(props1,2,mean)
p.qu <- t(apply(props1,2,quantile,prob=c(0.025,0.975)))

##############
# Scenario 2 #
##############

# Set covariates: Mixed conceptions and treatment 
covar2 <-  cbind(1, 0,1, 1, mean(attach), mean(pint), mean(pi_afd), mean(pi_fdp),
                 mean(pi_greens), mean(pi_left), mean(pi_spd), mean(pi_cducsu), mean(age),
                 mean(edulo), mean(eduhi), mean(female), mean(ost), 0,1)

# Calculate predictons 
props2 <- S%*%t(covar2)

# Means and confidence intervals 
p.mean2 <- apply(props2,2,mean)
p.qu2 <- t(apply(props2,2,quantile,prob=c(0.025,0.975)))

##############
# Scenario 3 #
##############

# Set covariates: One-sided conception, no treatment 
covar3 <-  cbind(1, 0,0, 0, mean(attach), mean(pint), mean(pi_afd), mean(pi_fdp),
                mean(pi_greens), mean(pi_left), mean(pi_spd), mean(pi_cducsu), mean(age),
                mean(edulo), mean(eduhi), mean(female), mean(ost), 0,0)

# Calculate predictons 
props3 <- S%*%t(covar3)

# Means and confidence intervals 
p.mean3 <- apply(props3,2,mean)
p.qu3 <- t(apply(props3,2,quantile,prob=c(0.025,0.975)))

##############
# Scenario 4 #
##############

# Set covariates: One-sided conception and treatment 
covar4 <-  cbind(1, 0,0, 1, mean(attach), mean(pint),mean(pi_afd), mean(pi_fdp),
                 mean(pi_greens), mean(pi_left), mean(pi_spd), mean(pi_cducsu), mean(age),
                 mean(edulo), mean(eduhi), mean(female), mean(ost), 0,0)

# Calculate predictons 
props4 <- S%*%t(covar4)

# Means and confidence intervals 
p.mean4 <- apply(props4,2,mean)
p.qu4 <- t(apply(props4,2,quantile,prob=c(0.025,0.975)))

##########
# PLOT   #
##########

# plot
exp.values1 <- c(props1, props2,props3, props4)
df1 <- data.frame(exp.values1)
df1$id <- c(rep("Control", 1000), rep("Treatment", 1000))
identity <- rep(1:2, times=1, each=2000)
df1 <- cbind(df1, identity)

labels <- c("2" = "One-sided conception","1" = "Mixed Conception")

#pdf("plot5.pdf")
ggplot(df1, aes(x=exp.values1, fill=id))+
  geom_density(alpha=0.4) +     
  guides(fill=guide_legend(title="")) + theme_bw() + scale_fill_manual(values=c("skyblue2", "dodgerblue4")) +                                  
  ylab("Density") + ggtitle("Effect of Framing Treatment on Support \nfor Helping Refugees", subtitle = "1=Agree, 5=Disagree")+
  xlab("Expected Value (1-5)") + xlim(2,3) + ylim(0,9) +facet_wrap(~identity, labeller=labeller(identity=labels)) +  
  theme(aspect.ratio = 1.2)
#dev.off()

#+ scale_fill_grey(start=0, end=0.6) +  

##################################
# Analysis with 2016 GLES data   #
##################################

rm(list=ls())

# Load GLES w1 data 
load("dataQ.Rdata")

############################
############################
# Models 1,2,3 of Table 4 ##
# + Figure 3              ##
############################
############################

modelhetGLES <- oglmx(w1_immi5~civ+ethn+edulo+eduhi+attach+w1_pint+
                female+age+ost, ~civ*ethn+attach+w1_pint+w1_pi_strength,
                link="logit", weights=weight1, constantMEAN=FALSE, constantSD=FALSE, data=dataQ)

summary(modelhetGLES)

modelhetGLESadd <- oglmx(w1_immi5~civ+ethn+edulo+eduhi+attach+w1_pint+
                   female+age+ost, ~civ*ethn+attach+w1_pint+w1_pi_strength+w1_uni*w1_trad,
                   link="logit", constantMEAN=FALSE,weights=weight1, constantSD=FALSE, data=dataQ)

summary(modelhetGLESadd)

modelhetGLESsimple <- oglmx(w1_immi5~civ+ethn+edulo+eduhi+attach+w1_pint+
                      female+age+ost,link="logit", constantMEAN=FALSE,
                      weights=weight1, constantSD=FALSE, data=dataQ)

summary(modelhetGLESsimple)

#####################
# Model fit         #
#####################

library("lmtest")

lrtest(modelhetGLESsimple,modelhetGLES)
lrtest(modelhetGLESsimple,modelhetGLESadd)

#############
# Figure 3  #
#############

coef <- coef(modelhetGLES)
coef <- as.matrix(coef)
vcov <- vcov.oglmx(modelhetGLES)
S <- mvrnorm(1000, coef, vcov)
Ss <- S[,10:15]
head(Ss)

##############
# Scenario 1 #
##############

# Set conception range
range <- seq(0,1, length.out=1000)

# Set covariates
values1 <- cbind(0.1,range, mean(attach), mean(w1_pint), mean(w1_pi_strength), 0.1*range)

# Check matrix
dim(values1)
head(values1)

# Multiply coefficients with values of independent variables 
props1 <- Ss%*%t(values1)

# Calculate mean and CI's
p.mean <- exp(apply(props1,2,mean))
p.qu <- exp(t(apply(props1,2,quantile,prob=c(0.025,0.975))))

###############
# Scenario 2  #
###############

# Set covariates
values2 <- cbind(0.9,range, mean(attach), mean(w1_pint), mean(w1_pi_strength), 0.9*range)

# Check matrix
dim(values2)
head(values2)

# Multiply with coefficients 
props2 <- Ss%*%t(values2)

# Calculate means and CI's
p.mean2 <- exp(apply(props2,2,mean))
p.qu2 <- exp(t(apply(props2,2,quantile,prob=c(0.025,0.975))))

##############
# Scenario 3 #
##############

# Set covariates
values3 <- cbind(range,0.1, mean(attach), mean(w1_pint), mean(w1_pi_strength), range*0.1)

# Check matrix
dim(values3)
head(values3)

# Multiply coefficients with values of independent variables 
props3 <- Ss%*%t(values3)

# Calculate mean and CI's
p.mean3 <- exp(apply(props3,2,mean))
p.qu3 <- exp(t(apply(props3,2,quantile,prob=c(0.025,0.975))))

###############
# Scenario 4  #
###############

# Set covariates
values4 <- cbind(range,0.9, mean(attach), mean(w1_pint), mean(w1_pi_strength), range*0.9)

# Check matrix
dim(values4)
head(values4)

# Multiply with coefficients 
props4 <- Ss%*%t(values4)

# Calculate means and CI's
p.mean4 <- exp(apply(props4,2,mean))
p.qu4 <- exp(t(apply(props4,2,quantile,prob=c(0.025,0.975))))

#############
# Make plot #
#############

# Append the different predictions and CIs
one <- list(p.mean, p.mean2,p.mean3,p.mean4)
one <- unlist(one)
two <- list(p.qu[,1], p.qu2[,1],p.qu3[,1], p.qu4[,1]) 
two <- unlist(two)
three <- list(p.qu[,2], p.qu2[,2],p.qu3[,2], p.qu4[,2])
three <- unlist(three)

# 4 times 0-1 range. 
x <- seq(0,1, length=1000)
s <- rep(x,4)

# 2 times 'low' and 'high' 
Conception <- rep(0:1,2, each=1000)
Conception <- ifelse((Conception==0),"Low", "High")

# Put data together 
identity <- rep(1:2, times=1, each=2000)
df.plot <- cbind.data.frame(one, two, three, s, Conception,identity)

labels <- c("1" = "Ethnocultural", "2" = "Civic")

#pdf("plot3_grey.pdf")
ggplot(df.plot, aes(x = s, y = one, ymax = three, ymin = two, group = Conception)) +
  geom_ribbon(alpha = 0.3, aes(fill = Conception))+ geom_line(aes())  + theme_bw()  + scale_fill_manual(name="Other \nConception", values=c("dodgerblue4","skyblue3")) +
  ylab("Standard Deviation of Residual") + xlab("Identity Conception") +
  ggtitle("Effect of National Identity Conception", subtitle="On Response Variability") + facet_wrap(~identity, labeller=labeller(identity=labels)) +  
  theme(aspect.ratio = 1.2)
#dev.off()

#####################################
# Analysis with GLES over-time data #
#####################################

rm(list=ls())

# load GLES over-time data
load("dataS.Rdata")

############################
############################
# Model 4 of Table 4      ##
############################
############################

#library(survey)
design1 <- svydesign(id=~1, weights=~weight2, data=dataS)
modeltest2 <- svyglm(sd2~civ*ethn+attach+w1_pint+w1_pi_strength+edulo+eduhi+
              w1_pi_afd+w1_pi_fdp+w1_pi_greens+w1_pi_left+w1_pi_spd+
              w1_pi_cducsu+female+age+ost, design=design1)

summary(modeltest2)  

# without weights
modeltest3 <- lm(sd2~civ*ethn+attach+w1_pint+w1_pi_strength+edulo+eduhi+
              w1_pi_afd+w1_pi_fdp+w1_pi_greens+w1_pi_left+w1_pi_spd+
              w1_pi_cducsu+female+age+ost, data=dataS)

summary(modeltest3) 

#####################################
#####################################
# Plot for the appendix            ##
# Marginal Effects for politically ##
# involved and not involved        ##
#####################################
#####################################

# Effect of ethnic conception over range of civic conception  

# Set range of civic
civic <- cbind(seq(0,1,0.1)) 

# coefficients 
coef <- coef(modeltest2)

# coefficients of interest
coef3 <- coef[3] # ethnic
coef18 <- coef[18] # civ*ethn

# Conditional interaction 
marg <- coef3 + (coef18*civic) 

# First, we need the variance of the estimated coefficients of interest:
Vpop <- vcov(modeltest2) # the Variance-Covariance matrix.
Vpop <- matrix(unlist(Vpop), ncol=18)
varpop <- diag(Vpop) # variance of all coefficients
var3 <- varpop[3] # variance ethnic 
var18 <- varpop[18] # variance interaction 

# we'll also need the covariance of the coefficients of interest:
cov <- Vpop[3,18] 

# Calculate variance of effect
var <- var3 + civic^2*var18 + 2*civic*cov 

# Now, we can calculate the standard errors:
se <- sqrt(var)

# With the standard errors, we can calculate the CIs as usual:
Ci.high <- marg + 1.96*se
Ci.low <-  marg - 1.96*se

#Prepare plot 
one <- marg
two <- Ci.high
three <- Ci.low
s <- seq(0,1, 0.1)
df.plot <- cbind.data.frame(one, two, three, s)

# Plot the conditional coefficients 

#pdf("s3.pdf")
ggplot(df.plot, aes(x = s, y = one, ymax = two, ymin = three)) +
  geom_ribbon(alpha = 0.3, fill="skyblue" ) + theme_bw() +
  geom_line()+
  geom_hline(yintercept=0, linetype=8, color="red") +
  ylab("Conditional coefficients") + xlab("Civic conception") +
  ggtitle("Effect of Ethnocultural Conception", subtitle = "Dependent variable = variance of immigration items (OLS)")
#dev.off()

# Set range of ethnic
ethnic <- cbind(seq(0,1,0.1)) 

# coefficients of interest
coef2 <- coef[2] # civic
coef18 <- coef[18] # civ*ethn

# Conditional interaction 
marg <- coef2 + (coef18*ethnic) 

# First, we need the variance of the estimated coefficients of interest:
Vpop <- vcov(modeltest2) # the Variance-Covariance matrix.
Vpop <- matrix(unlist(Vpop), ncol=18)
varpop <- diag(Vpop) # variance of all coefficients
var2 <- varpop[2] # variance civic 
var7 <- varpop[18] # variance interaction 

# we'll also need the covariance of the coefficients of interest:
cov <- Vpop[2,18] 

# Calculate variance of effect
var <- var2 + ethnic^2*var18 + 2*ethnic*cov 

# Now, we can calculate the standard errors:
se <- sqrt(var)

# With the standard errors, we can calculate the CIs as usual:
Ci.high <- marg + 1.96*se
Ci.low <-  marg - 1.96*se

# Prepare plot 
one <- marg
two <- Ci.high
three <- Ci.low
s <- seq(0,1, 0.1)
df.plot <- cbind.data.frame(one, two, three, s)

# Plot the conditional coefficients 

#pdf("sd4.pdf")
ggplot(df.plot, aes(x = s, y = one, ymax = two, ymin = three)) +
  geom_ribbon(alpha = 0.3, fill="skyblue" ) + theme_bw() +
  geom_line()+
  geom_hline(yintercept=0, linetype=8, color="red") +
  ylab("Conditional coefficients") + xlab("Ethnocultural conception") +
  ggtitle("Effect of Civic Conception", subtitle = "Dependent variable = standard deviation of immigration items (OLS)")
#dev.off()

# Other plots the same procedure without weights --> modeltest3

##################################
##################################
# Make plots for the appendix    #
# regarding political interest   #
##################################
##################################

##########################
# Political involvedness #
##########################

rm(list=ls())

load("data2.Rdata")

modelhet1 <- oglmx(ref3~civ+ethn+attach+pint+treat3+pi_afd+pi_fdp+pi_greens+pi_left+pi_spd+pi_cducsu+edulo+eduhi+age+female+ost,
                   ~civ*ethn*pint, data=data2, link="logit", weights=weight, constantMEAN=FALSE, constantSD=FALSE)

summary(modelhet1)

# subset data into politically involved and not involved 
data3 <- subset(data2, pint>0.5)
data4 <- subset(data2, pint<=0.5)

# among the politically involved 
modelhet2 <- oglmx(ref3~civ+ethn+attach+pint+treat3+pi_afd+pi_fdp+pi_greens+pi_left+pi_spd+pi_cducsu+edulo+eduhi+age+female+ost,
                   ~civ*ethn, data=data3, link="logit", weights=weight, constantMEAN=FALSE, constantSD=FALSE)

summary(modelhet2)

# among the politically not involved 
modelhet3 <- oglmx(ref3~civ+ethn+attach+pint+treat3+pi_afd+pi_fdp+pi_greens+pi_left+pi_spd+pi_cducsu+edulo+eduhi+age+female+ost,
                   ~civ*ethn, data=data4, link="logit", weights=weight, constantMEAN=FALSE, constantSD=FALSE)

summary(modelhet3)

################################################################### 
# CONDITIONAL COEFFICIENT PLOT FOR TRIPLE INTERACTION             #
###################################################################

###############################
# Political interest          #
###############################

# Range of ethnocultural
ethnic <- cbind(seq(0,1,0.1)) 

# coefficients 
coef <- coef(modelhet1)

############################
# Political interest low   #
############################

# coefficients of interest
coef1 <- coef[17] # civic
coef2 <- coef[18] # ethnic
coef3 <- coef[19] # polint
coef4 <- coef[20] # civ*ethn
coef5 <- coef[21] # civ*polint
coef6 <- coef[22] # ethn*polint
coef7 <- coef[23] # ethn*polint*civ

# Conditional interaction 
marg <- coef1 + (coef4*ethnic) + (coef5*0) + (coef7*ethnic*0)

# First, we need the variance of the estimated coefficients of interest:
Vpop <- vcov(modelhet1) # the Variance-Covariance matrix.
varpop <- diag(Vpop) # variance of all coefficients
var1 <- varpop[17] 
var2 <- varpop[18] 
var3 <- varpop[19] 
var4 <- varpop[20] 
var5 <- varpop[21] 
var6 <- varpop[22] 
var7 <- varpop[23] 

# we'll also need the covariance of the coefficients of interest:
cov14 <- Vpop[17,20] 
cov15 <- Vpop[17,21] 
cov17 <- Vpop[17,23] 
cov45 <- Vpop[20,21] 
cov47 <- Vpop[20,23] 
cov57 <- Vpop[21,23] 

# Calculate variance of effect
var <- var1 + ethnic^2*var4 + 0*var5 + ethnic^2*0*var7 + 
  2*ethnic*cov14  + 2*0*cov15 + 2*ethnic*0*cov17 + 2*ethnic*0*cov45 +
  2*0*ethnic^2*cov47 + 2*ethnic*0*cov57

# Now, we can calculate the standard errors:
se <- sqrt(var)

# With the standard errors, we can calculate the CIs as usual:
Ci.high <- marg + 1.96*se
Ci.low <-  marg - 1.96*se

# Prepare the plot 
one <- marg
two <- Ci.high
three <- Ci.low
s <- seq(0,1, 0.1)
df.plot <- cbind.data.frame(one, two, three, s)

#pdf("c3.pdf")
ggplot(df.plot, aes(x = s, y = one, ymax = two, ymin = three)) +
  geom_ribbon(alpha = 0.3, fill="skyblue" ) + theme_bw() +
  geom_line()+
  scale_y_continuous(limits = c(-3, 5)) +
  scale_x_continuous(limits = c(0, 1)) +
  geom_hline(yintercept=0, linetype=8, color="red") +
  ylab("Conditional coefficients") + xlab("Ethnocultural conception") +
  ggtitle("Effect of Civic Conception", subtitle= "Political interest = 0")
#dev.off()

############################
# Political interest high  #
############################

# coefficients 
coef <- coef(modelhet1)

# coefficients of interest
coef1 <- coef[17] # civic
coef2 <- coef[18] # ethnic
coef3 <- coef[19] # polint
coef4 <- coef[20] # civ*ethn
coef5 <- coef[21] # civ*polint
coef6 <- coef[22] # ethn*polint
coef7 <- coef[23] # ethn*polint*civ

# Conditional interaction 
marg <- coef1 + (coef4*ethnic) + (coef5*1) + (coef7*ethnic*1)

# First, we need the variance of the estimated coefficients of interest:
Vpop <- vcov(modelhet1) # the Variance-Covariance matrix.
varpop <- diag(Vpop) # variance of all coefficients
var1 <- varpop[17] 
var2 <- varpop[18] 
var3 <- varpop[19] 
var4 <- varpop[20] 
var5 <- varpop[21] 
var6 <- varpop[22] 
var7 <- varpop[23] 

# we'll also need the covariance of the coefficients of interest:
cov14 <- Vpop[17,20] 
cov15 <- Vpop[17,21] 
cov17 <- Vpop[17,23] 
cov45 <- Vpop[20,21]
cov47 <- Vpop[20,23] 
cov57 <- Vpop[21,23] 

# Calculate variance of effect
var <- var1 + ethnic^2*var4 + 1*var5 + ethnic^2*1*var7 + 
  2*ethnic*cov14  + 2*1*cov15 + 2*ethnic*1*cov17 + 2*ethnic*1*cov45 +
  2*1*ethnic^2*cov47 + 2*ethnic*1*cov57

# Now, we can calculate the standard errors:
se <- sqrt(var)

# With the standard errors, we can calculate the CIs as usual:
Ci.high <- marg + 1.96*se
Ci.low <-  marg - 1.96*se

one <- marg
two <- Ci.high
three <- Ci.low
s <- seq(0,1, 0.1)
df.plot <- cbind.data.frame(one, two, three, s)

#pdf("c4.pdf")
ggplot(df.plot, aes(x = s, y = one, ymax = two, ymin = three)) +
  geom_ribbon(alpha = 0.3, fill="skyblue" ) + theme_bw() +
  geom_line()+
  geom_hline(yintercept=0, linetype=8, color="red") +
  ylab("Conditional coefficients") + xlab("Ethnocultural conception") +
  ggtitle("Effect of Civic Conception", subtitle="Political interest = 1")
#dev.off()

# Repeat same for ethnic coefficient to produce plots showing marginal effect of
# ethno-cultural conception. 

#################################
# END OF REPLICATION FILE       #
# please contact Emmy Lindstam  #
# emmy.lindstam@gmail.com       # 
# for queries about additional  #
# analyses or raw data          #
#################################


